Do nonlinear waves in random media follow nonlinear diffusion equations? 
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Probably yes, since we find a striking similarity in the spatio-temporal evolution of nonlinear diffusion equations and wave packet 
spreading in generic nonlinear disordered lattices, including self-similarity and scaling. 
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1. Introduction 

The combined impact of disorder and nonlinearity strongly 
affects the transport properties of many physical systems lead- 
ing to complex behavior contrary to their separate linear coun- 
terparts. Application has great range; particularly relevant are 
nonlinear effects in cold atoms HI 12], superconductors [3], and 
optical lattices 0, [1, If!] . Still experimental probing both dis- 
ordered and nonlinear media remains limited due to reachable 
time or size scales. 

Significant achievements towards understanding the inter- 
play of disorder and nonlinearity have been made in recent 
theoretical and numerical studies. A highly challenging prob- 
lem was the dynamics of compact wave packets expanding in 
a disordered potential, in the presence of nonlinearity. The 
majority of studies focused on two paradigmatic models - 
the discrete nonlinear Schrodinger (DNLS) and the Klein- 
Gordon (KG) equations - revealing both destruction of an ini- 
tial packet localization and its resulting subdiffusive spreading, 
however with debate regarding the asymptotic spreading behav- 
iors Otis |9l lift ll llll2lll3ll . Hypotheses of an ultimate slow- 
down MM or eventual blockage of spreading lfl5ll have been 
recently challenged with evidence in 01611 . which reported a fi- 
nite probability of unlimited packet expansion, even for small 
nonlinearities. A qualitative theory of the nonlinear wave evo- 
lution in disordered media is based on the random phase ansatz 
[11], derives power-law dependencies of the diffusion coeffi- 
cient on the local densities, and predicts several distinct regimes 
of subdiffusion that match numerics convinc- 
ingly. Closely tied to these phenomena is thermal conductivity 
in a disordered quartic KG chain, analyzed in lfl7ll . 

Similar power-law dependencies of the diffusion coefficient 
on the local density have been extensively studied in the con- 
text of the Nonlinear Diffusion Equation (NDE). The NDE 
universally describes a diverse range of different phenomena, 
such as heat transfer, fluid flow or diffusion. It applies to gas 



flow through porous media 118, 19], groundwater infiltration 
lEol I2II1 . or heat transfer in plasma lEill . As a key trait, the 
NDE admits self-similar solutions (also known as the source- 
type solution, ZKB solution or Barenblatt-Pattle solution). It 
describes the diffusion from a compact initial spot and was first 
studied by Zel'dovich, Kompaneets, and Barenblatt I 23L I24I1 . 

The connection between nonlinear disordered spatially wave 
equations and NDE was conjectured recently and remains an 
open terrain lfl4l EM l26ll . A particularly challenging question 
is whether the NDE self-similar solution is an asymptotic time 
limit for the wave packet spreading in nonlinear disordered ar- 
rays. If yes, this will support the expectations that compact 
wave packets spread indefinitely, without re-entering Anderson 
localization. In this paper, we demonstrate that the NDE cap- 
tures essential features of energy /norm diffusion in nonlinear 
disordered lattices. At present, we still lack a rigorous deriva- 
tion of the NDE from the Hamiltonian equations for nonlinear 
disordered chains. Here we show that at sufficiently large time 
the properties of the NDE self-similar solution reasonably ap- 
proximate those of the energy/norm density distribution of non- 
linear waves; manifesting in similar asymptotical behaviors of 
statistical measures (such as distribution moments and kurto- 
sis), and in the overall scaling of the density profiles. To sub- 
stantiate our conclusions, we perform simulations of a modified 
NDE and compare the results against the spatio-temporal evo- 
lution of nonlinear disordered media 1 12l 131. 



2. Theoretical predictions 

2.1. Basic nonlinear disordered models 

The spreading of wave packets has been extensively stud- 
ied within the framework of KG/DNLS arrays. Particularly, 
the DNLS describes the wave dynamics in various experimen- 
tal contexts, from optical wave-guides ^ 01 to Bose-Einstein 
condensates B27I1 . The DNLS equation approximates well the 
KG one under appropriate conditions of small norm densities 
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1 28|, |29J . The KG has the advantage of faster integration at the 



same level of accuracy. 
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The DNLS chain is described by the equations of motion 

iifri = etfri + P\4ii\ 2 if/ 1 -iffi + \ - if/i-u (1) 

where e/ is the potential strength on the site /, chosen uniformly 
from an uncorrected random distribution [-W/2, W/2] param- 
eterized by the disorder strength W. 
The KG lattice is determined by 
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where m/ and pi are, respectively, the generalized coordi- 
nate/momentum on the site I with an energy density E\. The 
disordered potential strengths e; are chosen uniformly from the 
random distribution [1/2, 3/2]. The total energy E = Y,iEi acts 
as the nonlinear control parameter, analogous to ft in DNLS (see 
e.g. fiH). 

Both models conserve the total energy, the DNLS addition- 
ally conserves the total norm S — \if/i\ 2 . The approximate 
mapping from the KG to the DNLS is (3S w 3WE. We restrict 
analytics to the DNLS model, since the results are transferable 
to the KG one. 

2.1.1. Spreading predictions 

In order to quantitatively characterize the wave-packet 
spreading in Eqs. (1112b and compare the outcome to the NDE 
model, we track the probability at the l-th site, Pi = n\ = \if/[\ 2 , 
where n\ is the norm density distribution. Note that the analog 
of ni in the KG is the normalized energy density distribution 
E\. We then track a normalized probability density distribution, 
Zi = nil J]k n k- In order to probe the spreading, we compute the 
time-dependent moments m n = 2; zi(l - zf, where z — 2/ / zi 
gives the distribution center. 

We further use as an additional dynamical measure the kur- 
tosis lioll . defined as y(t) — m^if) / 'm£(i) - 3. Kurtosis is an 
indicator of the overall shape of the probability distribution pro- 
file - in particular, as a deviation measure from the normal pro- 
file. Large values correspond to profiles with sharp peaks and 
long extending tails. Low values are obtained for profiles with 
rounded/flattened peaks and steeper tails. As an example, the 
Laplace distribution has y = 3, while a compact uniform distri- 
bution has y = —1.2, 

The time dependence of the second moment ;«2 of the above 
distributions zi was previously derived and studied in 
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12l 1 1 311 . Different regimes of energy/norm subdiffusion were 
observed and explained. Generally, ma follows a power-law f 
with a < 1. Here we briefly recall the key arguments. 

In the linear limit Eqs. (1112b reduce to the same eigenvalue 
problem [l^l • We can thus determine the normalized eigen- 
vectors A v j and the eigenvalues A v . With if// = 2Zv^v,l<Pn Eq.© 
reads in an eigenstate basis as 



V,Vl,V2,V3'rvi T» 2 rV3 J 



(3) 



where I v , VilV%V i - TiiA Vt iA vu iA v%i iA Vi> i are overlap integrals and 
<p v determine the complex time-dependent amplitudes of the 
eigenstates. 



In Ultl the incoherent "heating" of cold exterior by the 
packet has been established as the most probable mechanism 
of spreading. Following this analysis, the packet modes <p v {t) 
should evolve chaotically with a continuous frequency spec- 
trum. In particular chaotic dynamics of phases is expected 
to destroy localization. The degree of chaos is linked to the 
number of resonances, whose probability becomes an essen- 
tial measure for the spreading. Previous studies 113111 indicate 
that the probability of a packet eigenstate to be resonant is 
Kifin) — 1 - e~ c/3 ", with C being a constant dependent on the 
strength of disorder. The heating of an exterior mode close to 
the edge of the wave packet with norm density n would then 
follow 

% = W M + /3n 3/2 'R(j3n)f(t) (4) 

with delta-correlated (or reasonably short-time correlated) 
noise /(?), and lead to UJ ~ y8 2 n 3 ('R(/?n)) 2 f. The momentary 
diffusion rate follows as D ~ /3 2 n 2 CR(J3n)) 2 . 

With m 2 ~ n~ 2 = Dt one arrives at 1 /n 2 ~ {3(1 - e~ ci3n )t 112 . 
Depending on the product C/3n being larger or smaller than one, 
the packet has two regimes of subdiffusion (and a dynamical 
crossover between them): ;«2 ~ fit 1 / 2 (strong chaos) and ni2 ~ 
ff'h 113 (asymptotic weak chaos) iHllSlIIIllQJ]. 

The validity of the assumption of incoherent phases and of 
Eq.© was established through numerical studies for the first 
time by Michaely and Fishman 113211 . moving the above conjec- 
ture based theories onto solid grounds. 

2.2. Nonlinear diffusion equation 

The assumption of chaos and random phases, Eq.©, the den- 
sity dependent diffusion coefficient and the resulting subdiffu- 
sion strongly suggest an analogy to nonlinear diffusion equa- 
tions (see e.g. Ref.[33]). We first consider here the much stud- 
ied NDE with a power-law dependence of the diffusion coef- 
ficient on the density. The NDE in the one-dimensional case 
reads 

d t 9 = d x (Kf a d x f f ) (5) 

Here P = f(x, t) is the concentration of the diffusing species 
(which may be related to the energy/norm density), a > and 
k are some constants. Hereafter, we set k — 1 without loss of 
generality. 

Let us discuss the scaling properties of Eq.©. Given a solu- 
tion fix, t) it follows that 



P(x, f) = s p P(s x x, s,t) 



(6) 



where s p , s x , s t are some scaling factors. From normalization it 
follows that s p = s x . Inserting © into (© we find 

s, = C 2 ■ (7) 
After some algebra for the moment m n {t) we conclude 

\ >?/(a+2) 



mJt ) 



(8) 



which corresponds to a subdiffusive process. 
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Figure 1: The NDE parameter a versus the energy density E{t) ~ E(0) (ni2(0) '' 2 as obtained from numerical simulations of the KG chain. Horizontal lines guide 
the eye at values a = 2 and a = 4. Left panel: W = 4, L = 21, and E e {0.01, 0.015, 0.02, 0.04, 0.06, 0.1, 0.75) (varying from the top to bottom and from weak chaos 
to strong chaos and self trapping). Right panel: W = 2 and E = 0.1 (curve 1, strong chaos), W = 4 and E = 0.2 (curve 2, crossover from the strong to weak chaos), 
W = 4 and E = 0.01 (curve 3, weak chaos), W = 6 and E = 0.05 (curve 4, weak chaos). Inset: Dependence of a on time for the same data as in the right panel. The 
dot-dashed and dashed lines correspond to the values a = 1 /3 and a = 1/2 respectively. 



Eq.© has various self-similar solutions 11341 l35l Kill which 
depend on the characteristics of the evolving state. For compact 
wave packets the self-similar state realized in the asymptotic 
limit of large time is of the form 11231 l24l 13611 
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where A is a normalization constant and |x| < xq. For |x| > xq 
the density P strictly vanishes. The edge position xo depends 
on time as 

x = [2A(2/a + l)f 2/(2+a) ]' . (10) 

The linear stability of Eq.© was demonstrated in |37, 38]. 

Using a change of variables y = x/xo we obtain for the den- 
sity P(y, t) with Pdy = Pdx: 



p(y, t ) = A I/fl+1/2 V4/a + 2(1 - y 2 ) l/a 
Since P does not depend on time, it follows 



(11) 



1,(0 = x\m n , m n = f yipdy . (12) 



In agreement with (0 this yields e.g. mi 
The moments of solution ((9) 
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where S(x,y) is the Euler Beta Function. Using Eq.(fT3l we 
derive for the kurtosis 



joo = lim y = 3 



-3, 



(14) 



where is Legendre's Gamma function. As can be seen, the 
kurtosis of the self-similar solution does not depend on time. 
For the values a — 4 and a — 2 we obtain kurtosis values y = 
-1.091 and y = -1.00 respectively. Additionally, y = -1.20 
in the limit of a — » oo, which corresponds to a flat uniform 
distribution. 

Few remarks are in place in order to proceed. First, the 
spatial discretization of the NDE (0 by introducing discrete 
Laplacians with nearest neighbor differences does not modify 
the properties of the asymptotic states [33]. However, the over- 
lap integrals in Q decay exponentially in space. Moreover, the 
diffusion coefficient D for spreading wave packets is in general 
different from a pure power of the density (see (|4} and below). 
It takes a power function only in the asymptotic regime of weak 
chaos, and the potential long lasting intermediate strong chaos 
regime. Therefore we will generalize and adapt the above NDE 
in the next subsection. 

2.3. Modified nonlinear diffusion equation 
Firstly, we rewrite Eq.(|5) as 
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Since density leakage into nei ghbo rin g si tes is directly related 
to resonance probability Ir9l ll0l.ll 1 . 12, 13], we introduce it into 
the RHS ofEq.Oas 



Randomness in disordered systems exponentially localizes the 
normal modes, so mode-mode coupling in the nonlinear overlap 
integral has an exponential dependence in distance as well. We 



therefore introduce an exponentially decaying interaction along 
a discrete chain. Using a finite central difference, we arrive at 
the modified nonlinear diffusion equation (MNDE) 



p -mlx (p _ 2F + F , ) 
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(16) 



In the above, we treat C as a free parameter, and equal x m our 
numerics to localization lengths of the disordered lattice mod- 
els. We fix a — 2 for the MNDE, and expect to observe similar 
regimes, such that for Pi ~ 1, we have weak chaos for C <sc 1 
and strong chaos in the opposite limit of C » 1. The above 
MNDE is therefore expected to account for the resonance prob- 
abilities between normal modes, and the exponentially decay- 
ing interaction between them. 

2.4. Numerical simulations of MNDE 

We integrate Eq.([T6l> with a fourth-order Runge-Kutta 



scheme 03911 . for a number of values for the free parameter C. 
We start with an initially compact distribution of width L = 41 
and density P\h<l/2 = 1 (hereafter referred to as brick distri- 
bution) and !P|;|>l/2 = 0. The integrations were carried out to 
t « 10 6 using a time step of 0.4, all the while conserving norm 
to better than 1(T 12 . 

From the second moment, m2(t), we compute the derivative 



a(t) 



dlog 10 OT 2 
dlog 10 f 



(17) 



and plot the result in the upper panel of FigfJ] We find that 
the MNDE reproduces weak chaos (a = 1/3 for C < 1) and 
intermediate strong chaos {a - 1/2 for C » 1). 

We further plot in Figf2]the time-dependent kurtosis (lower 
panel) and the density distributions at t — 10 6 (inset) for a few 
representative values of C. Note that those states that are in the 
weak chaos regime (C < 1) show a tendency toward an asymp- 
totic yoo » -1.091 in agreement with the NDE, but not reach- 
ing it fully in our simulations. Those states in the intermediate 
strong chaos regime (C > 100) exhibit long lasting saturation at 
7oo ~ -1 in agreement with the NDE. We also observe a growth 
of the kurtosis into positive values for weak chaos, followed by 
a drop that decays to - 1 .09 1 . 

3. Wave packet spreading in nonlinear disordered lattices 

Let us discuss first the details of computations. For both 
models of Eqs. (lll21 > we consider initial brick distributions of 
width L with nonzero internal energy density (or norm density 
for the DNLS) and zero outside this interval. In contrast to the 
NDE and MNDE, these lattice systems, in addition to local den- 
sities, are also characterized by local phases. Initially the phase 
at each site is set randomly. Equations ([TJ, (f2]i are evolved us- 
ing SABA-class split-step symplectic integration schemes [40], 
with time-steps of 10" 2 - 10 _I . Energy conservation is within 
a relative tolerance of less than 0.1%. We perform ensemble 
averaging over 10 3 realizations of the onsite disorder. 

With m 2 ~ t 2/(2+a) of the NDE self-similar solution, Eq.(fT3l. 
and 1112 ~ f for KG/DNLS models, the NDE parameter a is 
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Figure 2: Numerical results for the MNDE, Eq. dl6t with x — 6 and a = 2 and 
initial conditions given in the text. Upper: The exponent a(t) for the MNDE 
with the following values of C: 0.1 (dashed black), 1.0 (solid black), 10 (dashed 
dark grey), 25 (solid dark grey), 50 (dashed light grey), and 1000 (solid light 
grey). Lower: The kurtosis y(i) for the same simulations. Inset: The density 
profiles of the above simulations at t = 10 s . 



related to the exponent a as a = 2(1 - a) I a. This allows a 
monitoring of a as the energy density changes. We expect then 
a — 2 and a — 4 respectively for the strong and weak chaos 
regimes, as well as a shift from a — 2 to a — 4 associated with 
the crossover between the two regimes. 

We validate these predictions using our numerical data that 
correspond to the different regimes of spreading [12]. We 
plot the NDE parameter a from the numerically obtained a 
(see inset in FigQ]), assuming that the energy density E{t) ~ 
E(0) (/tt2(£)) 2 in FigQ] As predicted, a reaches the asymp- 
totic value a = 4 for weak chaos. Our numerical results also 
show that once a reaches its asymptotic value a = 4, it does 
not increase further in time, even for quite small energy den- 
sities. That is a clear indication for the absence of speculated 
slowing down dynamics [Mfl. For strong chaos a temporarily 
saturates around a — 2, keeps this value only within some inter- 
val of energy densities, and finally crosses over into the interval 
2 < a < 4 with a clear tendency to reach the weak chaos value 
a — 4 at larger times. 

The resulting kurtosis evolution, (y(f)) is presented in Figf3] 
For the initial wave packet (y(0)) = -1.2. The kurtosis first dis- 
plays a transient increase to positive values. This is very similar 
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Figure 3: Semi-log plot of average over 10 3 realizations kurtosis (y) versus 
time for KG (see main part) and DNLS (see inset) models with parameters 
W = 4, L = 21. Numbers correspond to different energy densities (KG), or, 
nonlinearity strengths (DNLS): E = 0.01 or p = 0.04 (curve 1), E = 0.02 or 
P = 0.08 (curve 2), E = 0.04 or p = 0.18 (curve 3), E = 0.08 or p = 0.36 
(curve 4), E = 0.2 or p = 0.72 (curve 5). The dashed lines correspond to the 
<y> = -1.0. 



to the MNDE results and is due to exponential localization of 
the initial state in normal mode space. At larger times (y(t)) 
displays a decrease in time, approaching the self-similar behav- 
iors in density distributions with (y) « — 1 (recall that the NDE 
self-similar solution gives us y — -1 for a — 2 and y = -1.091 
for a = 4). 
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Figure 4: Main: The log of the normalized energy density distribution (log 10 zi) 
at three different times (from the top to bottom t a 10 4 , t a 10 7 , t a 10 s ). The 
initial parameters are E = 0.2, W = 4 and L = 21. Symbols correspond 
to the average over 10 3 disorder realizations, and solid lines correspond to an 
additional smoothing. Inset: Rescaled distributions (see text). 

The evolution of the averaged energy density profiles (E) in 
the course of spreading is illustrated in Figj4] The peaked ini- 
tial distribution profiles transform into more fiat ones as time 
evolves. The most striking result is obtained by rescaling the 
profiles in Fig |4] according to the scaling laws of the NDE. We 
estimate the value of xq (see (fTUl l) as the position where the 
profiles in Figj4] reach 10~ 4 . We then plot the rescaled den- 



sities according to (TPTT i in the inset of Fig|4] We observe very 
good scaling behavior. This is the strongest argument to support 
the applicability of NDE and MNDE to the spreading of wave 
packets in nonlinear disordered systems. It also strongly sup- 
ports that the spreading process follows the predicted asymp- 
totics and does not slow down or even halt. 

4. Conclusion 

We scrutinized the suggested connection between the tempo- 
ral evolution of self-similar solutions of the NDE and MNDE 
on one side, and the asymptotic dynamics of the energy/norm 
distribution within nonlinear disordered media on the other, 
and found a remarkable correspondence. In order to describe 
the expansion of an initial distribution with time, we used two 
key quantities: (i) the second moment m2(t), which shows how 
the squared width of the distribution grows; (ii) the kurtosis 
y(t), which indicates how the shape of the distribution profile 
changes. 

As a first test we compared the exponents characterizing the 
subdiffusion for the second moments of the energy/norm den- 
sity distributions. In KG/DNLS models, m2(t) ~ f , and for 
the NDE self-similar solution, TO2(?) ~ ? 2 / (2+fl) , the exact iden- 
tity giving a = 2/(2 + a). We found that the wave packets in 
nonlinear disordered chains converge towards the self-similar 
behavior at large times. The numerical results show a good 
correspondence to the NDE-based analytics in a wide range of 
parameters. 

Second, in disordered lattices the energy/norm distributions 
have exponentially decaying tails at variance to the steep-edged 
NDE self-similar solution. Such difference of energy/norm dis- 
tribution profiles has no effect on ni2(t) at large times. However, 
it leads to differences between the NDE and the KG/DNLS dy- 
namics at intermediate times, e.g. seen in the temporal behavior 
of the kurtosis (y(f)). In order to study the possible impact of 
various initial density profiles, we also computed the time evo- 
lution of the NDE for initial probability distributions with expo- 
nential tails. In all simulations with NDE parameter a < 6, the 
kurtosis asymptotically reached the expected value y^ being a 
function of the parameter a only. 

To bridge a possible gap between the NDE and the disor- 
dered nonlinear lattices, we introduced a modified MNDE. To 
account for the interaction between localized Anderson lat- 
tice modes we implemented exponentially decaying coupling 
in the MNDE, and also incorporated the resonance probabil- 
ity of normal modes into a modified nonlinear diffusion coef- 
ficient. Then we indeed observe the dynamical behavior that 
reproduces the spatio-temporal evolution of nonlinear disor- 
dered chains, namely the weak and the strong chaos regimes 
of spreading, and he correct temporal evolution of the kurtosis, 
and the distribution profiles. Most importantly we observe the 
precise scaling behavior of the asymptotic NDE solutions in the 
case of nonlinear wave packets. 

Let us summarize. There is a lack of knowledge on the sta- 
tistical properties of chaotic dynamics generated by nonlinear 
coupling. We are still far from a rigorous derivation of the NDE 
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when starting with the equations of motion for nonlinear disor- 
dered chains. Nevertheless, the theory for initial energy/norm 
spreading in KG/DNLS chains which is based on a Langevin 
dynamics approximation has earlier been confirmed by exhaus- 
tive numerical studies. Of course, there is difference between 
KG/DNLS nonlinear disordered models and the NDE. Despite 
this, numerical results confirm that at sufficiently large time 
the NDE self-similar solution approximates remarkably well 
the spreading properties of energy/norm density distributions 
in terms of the second moment, the kurtosis, and the scaling 
features. Therefore, the NDE as a simple analytical tool is ex- 
tremely useful for studying the initial excitation spreading in 
nonlinear disordered media at asymptotically large times. Ad- 
ditionally, the NDE analog with long-range exponentially de- 
caying coupling shows an even deeper correspondence between 
generic nonlinear disordered models and the NDE and therefore 
might prove to be an insightful model for the future analysis of 
spreading in nonlinear disordered systems. 

We thank N. Li, M.V. Ivanchenko, and G. Tsironis for in- 
sightful discussions. 
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